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ABSTRACT. In this paper, we examined the effects of interactions in a predator- 
prey manner and harvesting at constant fishing efforts of Nile perch (Lates 
niloticus), Nile tilapia (Oreochromis niloticus) and Small pelagic silver (Ras- 
trineobola argentea) fish species of the lake Victoria fishery, through mathe- 
matical modeling using Lotka-Volterra equation, whereby all three fish species 
are subjected to harvesting. The model is characterized by the system of first 
order non-linear ordinary differential equations. 

All eight equilibrium points of the model were identified, the local stability of 
the co-existence equilibrium point was discussed using Routh-Hurwitz criteria 
and its global stability was analyzed using suitable Lyapunov function. Fur- 
ther, analytic solutions of the model coincided with the computed numerical 
solutions using fourth order Runge-Kutta method. The study revealed that as 
the model parameters became small the equilibrium stock level biomasses of 
fish species increased. 


1. INTRODUCTION 


Lake Victoria is the second largest fresh-water lake in the world and the largest 
in Africa. The lake is shared by Kenya (6%), Uganda (43%) and Tanzania (51%). 
According to [7], lake Victoria fishery is currently dominated by three fish species 
which are Nile perch, Nile tilapia and small pelagic silver fish. The lake is an im- 
portant source of food, employment and earnings for the riparian states through 
fishery for millions of people. 


The study aims at investigating the effects of harvesting at constant fishing ef- 
forts and interactions in a predator-prey manner to the equilibrium stock level 
biomasses of Nile perch, Nile tilapia and small pelagic silver fish species in lake 
Victoria. The Nile perch are predator to Nile tilapia and vice versa (depending 
on their sizes, that is, adult Nile perch eats the young Nile tilapia and adult Nile 
tilapia eats the young Nile perch) while small pelagic silver fishes are prey to both 
Nile perch and Nile tilapia. 

Quantitative and qualitative understanding of the interactions of these fish species 
is crucial for the management of fisheries in lake Victoria. Harvesting has generally 
a strong impact on the population dynamics of harvested species. The severity of 
this impact depends on the nature of the applied harvesting strategy which, in turn, 
2010 Mathematics Subject Classification. 92-08, 65P99, 92B05, 37N25. 
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may range from the rapid depletion to the complete preservation of a population|[10]. 


The problem of interspecies interactions was considered by [4] for two species 
obeying the law of logistic growth. [2] considered harvesting of a single species 
in an ecologically competing two fish population model.{1] and [11] studied the 
dynamics of two-species fishery by combining harvesting. [3] studied constant rate 
of harvesting in a predator-prey system to allow simultaneous harvesting of both 
species. They showed how to approximate the region of asymptotic stability in 
biological terms, in the initial states which lead to coexistence of the two species 
and their global dynamics by efficient computer simulation. [9] studied the effect 
of Nile perch predation to Nile tilapia and harvesting on fisheries dynamics in Lake 
Victoria, the study ignored the third species, the small pelagic silver fish, a prey 
to the first two species which has a significant contribution to the dynamics of 
species.To the authors’ knowledge,the nature of this study has not yet been done. 


2. MATHEMATICAL MODEL 


2.1. Assumptions of the model. The model rely on the following assumptions: 


e The fish species can grow independently in a lake and their population sizes 
are bounded 

e Adult Nile perch are predator to both young Nile tilapia and Small pelagic 
silver fish 

e Adult Nile tilapia are predator to both young Nile perch and Small pelagic 
silver fish 

e Small pelagic silver fishes are prey to both Nile tilapia and Nile perch 

e All three fish species have ecological interactions in a lake 

e All three fish species are harvested at a rate proportional to the size of their 
population 

e Fishing effort for all fish species is kept constant 


2.2. Definitions of variables and parameters of the model. The fol- 
lowing are the definitions of variables and parameters used in developing 
the model: 
rı: intrinsic growth rate of Nile perch 
rg: intrinsic growth rate of Nile tilapia 
r3: intrinsic growth rate of small pelagic silver fish 
: stock biomass of Nile perch 
: stock biomass of Nile tilapia 
stock biomass of small pelagic silver fish 
: predation rate of Nile tilapia to Nile perch 
: predation rate of Nile perch to Nile tilapia 
: predation rate of Nile perch to small pelagic silver fish 
: predation rate of Nile tilapia to small pelagic silver fish 
E: fishing effort for Nile perch 
FE: fishing effort for Nile tilapia 
Ez: fishing effort for small pelagic silver fish 
qı: catchability coefficient of the Nile perch 
q2: catchability coefficient of Nile tilapia 
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q3: catchability coefficient of small pelagic silver fish 
ky: carrying capacity of Nile perch 

Kə: carrying capacity of Nile tilapia 

Ks: carrying capacity of small pelagic silver fish 


2.3. Model equations. The model equations for the three fish species is 
the set of non-linear ordinary differential equations: 


dx 
q Ele at Py +72) 
d 
= = y (co + px — azy + 42) 
t 
dz 
q 72 (8 — WE Vy — asz) 


where a; = x > 0 for i = 1,2,3 ,¢c; = rj —qj;E; > 0 for j = 1,2,3 and 
(a—B)=p>0 


3. EXISTENCE OF EQUILIBRIUM POINTS 
The system under investigation has eight possible equilibrium points ob- 


d d 
tained by setting ai = a = < 


(i) Fy (x*, y*, 2*) = (0,0,0) (absence of all fish species) 
(ii) Ea (a*,y*, 2*) = (o 0, =) (absence of Nile perch and Nile tilapia) 
a3 


= 0.These includes: 


(iii) Es (a*,y*, 2") = (2, 0, 0) (absence of Nile tilapia and small pelagic 
ay 
ilver fish) 
(iv) E4 (a*,y*,2*) = (o 2 0) (absence of Nile perch and small pelagic 
a2 


n 


silver fish) 
x oe —Cop + A2C, Q1C2 + Cp 
(v) Es (a*,y", 27) = 2 179 
pP + a,a2 pP + a,ag 
silver fish). The equilibrium point Es is positive if agc, > cep hold. 


(vi) Es (a*,y*,2*) = = as 2 ; ae ante il (absence of Nileti- 
Y“ + ay ag y+ a1 a3 

lapia).The equilibrium point Ee is positive if a1c3 > c1y hold. 
(vii) Ez (2*, y*, z*) = 0, Wc3 T aed a2Cc3 — cow 
Y? + a2a3° Y? + azaz 

The equilibrium point E7 is positive if agc3 > c2% hold. 





: 0) (absence of small pelagic 











(absence of Nileperch). 

















—pre3—a3 pc2—ywoe2 tyc3a2 +e1y)?+4a2403¢1 
a Z 2 2 
azp? +a% +a1a2a3+7°a2 
(väi) Bs ( yt |=| aamtapeteeter -neten 
p 3P 1 142a3 TY 2, 
zZ —(—Ypc2+yc1a2+Yc201 +YpcI -c3 p“ —c3a102) 
a3p?+aı p? +a1a2a3+7%?a2 











(the co-existence of all three fish species).The equiliribrium point Eg is 
positive provided the following inequalities holds: 
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Yczaz + Y? + azazcı > pivc3 + aspc2 + yWre2, 
a1 a3C2 + agcip + pyc3 + c27? + Ycsaı > ye, and 
ypco + esp” + c3a1a2 > ye1a2 + Wea + Ype 





3.1. Local stability analysis. The local stability of the interior equi- 
librium point (co-existence point Eg) was investigated using the Routh- 
Hurwitz’s criteria. That is, an equilibrium point is locally asymptotically 
stable if the characteristic equation of the Jacobian matrix evaluated at 
that point has all the coefficients being positive and that all of its roots 
have negative real parts[8]. 


J(Es) = 
cı — 2a,2* — py* + 2" —px* yr* 
py* C2 + px* — 2a2y* + pz" vy" 
aye —yz* c3 — yx" — py* — 2a3z* 
The characteristic equation of J (Eg) is 
(3. 4) Nga (by + bə + b3)? + (by be + bob3 + bibs) = (bi b2b3 + ba) =0 
where 
by = cy — 2a,x* — py“ + yz* = —aıx* < 0 since, 2* > 0 and 





cy — ayx* — py* + yz* =0 

bo = co + px* — 2a2y* + wz* = —a2y* < 0 since, y* > 0 and 
C2 + px* — agy* + %z* =0 

b3 = c3 — ya* — wy* — 2a3z* = —a3z* < 0 since, z* > 0 and 
c3 — yu* — wy* — a3z* =0 

by = wy*z* > 0 

Hence 

(by + b2 + bs) <0> — (by + b2 + bs) >0 

Similary, (by be + bob3 + by bs) >0 

And 











bibeb3 +b4 = —azagagx*y*z* + y2y*z* 
= y*2*(W? — ayaza32°*) 
Certainly, ajaza3x* > y? 
Hence, bi bob3 + b4 = y*2* (Y? — ayaza32*) < 0 = — (bi bab3 + b4) > 0 
The conditions for Routh-Hurwitz’s criteria are satisfied, that is, 
— (bı + bə + bs) (b1b2 + bob3 + bbs) > —(b1b2b3 + ba) 


Therefore, the co-existence equilibrium point Eg is locally asymptoti- 
cally stable. 


3.2. Global stability analysis. Global stability of the system was anal- 
ysed by considering suitable Lyapunov function [5]. 
Consider the following Lyapunov function candidate, 


(3. 5) V(a,y, 2) =h[e@—2*—-2* m(Ž)] Holy—y*—y* In) sl2 22" In(=)] 





where 11, lo, 13 > 0. 


(3. 6) 


(3. 7) 
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As per, [6], it is evident that the choosen Lyapunov function candi- 
date V (x,y,z) of ( 3. 5 ) satisfy the conditions that V(a*,y*,z*) = 0 and 
V(a,y,z) > 0 for all (x,y,z) 4 (a*, y*, 2*). Moreover, V(x, y, z) is radially 


unbounded. 


We are required to verify that 


0,l2 > 0 and l3 > 0. 





av 
dt 


< 0 for the suitable choice of l > 



























































E TE 10-2 
= Eg HOS ag Sg 
Let 
~ =F+G+H 
where, 
F = lı (x — z*)(cı — aix — py + yz) 
G = lly — y*)(c2 + px — azy + 4z) 
H =13(z — z*)(c3 — yz — Wy — a32) 
Considering, 
F = h(x- zr*)(c& — az — py + 72) 
= _ hle(a — az — py +yz) — x* (c1 — ax — py + 72)] 
= hlar- a,x? — pry +yrz— cr* + agzz* + pyx* — yzx*] 
Since, cı — a1£* — py* + yz* = 0, > c1 = ay 2* + py* — yz" 
Then 
F = h[e(az* + py* — y2*) — az’ — pry + yxz 
x*(a1£* + py* — yz") + ayaa* + pyx* — yzx*] 
= h[-ai(2* — 2ea* +2") — ply(x—2*) —y* (a — 2*)| 
Ho y[—2"*(a~ — 2*) 4+ 2(a — 2*)]] 
Further algebraic simplification gives, 
F = lı (x — 2") [~ai (z — 2*) — p(y — y*) + (2 — 2*)] 
Considering, 
G = l(y-— y*)(c2 + px — azy + oz) 
= 12 [y(co + px — azy + Yz) — y* (c2 + px — azy + ¥2)] 
= lə (coy + pry — a24? + pyz — cay* — pry* + aayy* — py“ z) 
Since, cz + px* — azy“ + pz* = 0, = co = —pa* + azy“ — pz" 
Then 
G = lely(—px* + azy" — pz") + pry — azy’ + pyz 
y"(—px" + aay” — p2") — pry” + azyy* — py" z] 
= ly[-a2(y—y*)? + p[-a*(y—y*) + 2(y—-y"*)] 
t pl-2"(y—y") + 2(y—y*)I] 
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Further algebraic simplification gives, 





(3. 8) G = lə(y — y*) [p(x — x“) — aa(y— y*) + Y(z — 2*)] 
Considering, 
H = lg(z— z*)(c3— ya — py — azz) 








= la [z(cg — yz — wy — agz) — z* (c3 — yz — Wy — a32)] 


l3 [zc3 yrz — wyz — azz? — c3z* + yar* + yz" +4 azzz“] 





Since, c3 — ya* — py* — a3z* = 0, > c3 = ya* + Wy* + a3z* 
Then 





H = Islz(ya* +Wy* + a3z*) — yaz — yz — azz? 
— z*ř(yx* + wy* + agz*) + yxz* +yz“ + a3zz“] 
= Is[—ag(z — 2*)? — y(x — 2*)(z — 2*) 











Further algebraic simplification gives 





(3. 9) H = I3(z— 2") [-7(@— 2") — p(y —y") — as(z — 2") 


Substituting (3. 7 ),(.3. 8 ) and ( 3. 9 ) into ( 3. 6 ) gives, 


dV 
(3. 10) es 1X (—a,.X —pY +yZ)+leY (pX —ag¥ +Z)+13Z(—yX -4Y —a3Z) 





where, X = (x — 2*), Y = (y—y*) and Z = (z — 2") 


Further simplification of ( 3. 10 ) gives the following: 





So SS [lia X? + bay? + lza3Z°] t loll L)XY H v(h 13)XZ + W(lg — l3)Y Z] 


If (X,Y,Z) = (0,0,0), that is, when [x = z*, y = y* and z = 2*] then 
dV 
= 0. 


dt 
And, if l = l2 = Is, then w <0, for all (2*, y*,2*) A (x,y,z) in R3. 
Therefore, the co-existence equilibrium point Eg is globally asymptoti- 
cally stable. 


4. NUMERICAL EXAMPLES 


Numerical examples and their graphical illustrations are summarized 
below in eight different cases. Fourth order Runge-Kutta integration algo- 
rithm was used for some cases to validate the qualitative analysis results. 


= a ma w 
A a 0o [=] 


species population 


tas 
bo 
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TABLE 1. Parameters for figure la and figure 1b of co-existence 


equilibrium point Eg 









































Parameters | Figure la | Figure 1b 
ay 0.0125 0.003 
ag 0.07 0.055 
a3 0.01 0.001 
p 0.04 0.02 
y 0.005 0.001 
Y 0.005 0.001 
Cl 0.95 0.91 
C2 0.50 0.70 
C3 0.40 0.21 
(a*,y*,2*) | (20,20,20) | (70,40,100) 
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FIGURE 1. Graphical representations of the parameters in Table 1 


4.1. Examples 1&2 with z(0) = 10, y(0) = 10 and z(0) = 10 for both 


cases. 


TABLE 2. Parameters for figure 2a and figure 2b of co-existence 


equilibrium point Eg 















































Parameters Figure 2a Figure 2b 
ay 0.0015 0.002 
az 0.005 0.0035 
a3 0.001 0.0005 
p 0.002 0.0002 
y 0.001 0.0005 
w 0.001 0.0005 
C1 0.75 0.22 
C2 0.60 0.22 
C3 0.90 0.225 
(x*,y*,z*) | (300,300,300) | (150,100,200) 
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FIGURE 2. Graphical representations of the parameters in Table 2 


4.2. Examples 3&4 with x(0) = 10, y(0) = 10 and z(0) = 10 for both 


cases. 


TABLE 3. Parameters for figure 3a and figure 3b of co-existence 


equilibrium point Eg 















































Parameters Figure 3a Figure 3b 
ay 0.0015 0.0008 
az 0.003 0.001 
a3 0.004 0.0005 
p 0.002 0.0003 
y 0.0035 0.0004 
w 0.001 0.0004 
Cy 0 0.42 
C2 0 0.18 
C3 0.85 0.78 
(x*,y*,z*) | (100,100,100) | (600,600,600) 
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FIGURE 3. Graphical representations of the parameters in Table 3 
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4.3. Examples 5&6 with (0) = 10, y(0) = 10 and z(0) = 10 for both 
cases. 


TABLE 4. Parameters for figure 4a and figure 4b of co-existence 
equilibrium point Es 






































Parameters Figure 4a Figure 4b 
ay 0.000013 0.00000024 
az 0.000015 0.00000022 
a3 0.000004 0.00000008 
p 0.000001 0.00000004 
y 0.000004 0.00000008 
Y 0.000004 0.00000008 
Cy 0.10 0.020 
C2 0.10 0.010 
C3 0.12 0.024 
(x*,y*,z*) | (10000,10000,10000) | (100000,100000,100000) 
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FIGURE 4. Graphical representations of the parameters in Table 4 


4.4. Examples 7&8 with x(0) = 50, y(0) = 10,z(0) = 10 and z(0) = 
10,y(0) = 10,z(0) = 10 respectively. 
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TABLE 5. Runge-Kutta numerical results for the parameters of figure 1b 














t x(t) y(t) z(t) 
0 10 10 10 
1 18.9543128777437318 14.0106457936006752 11.8888854688601402 
2  31.8253617861179486 19.1105969078411135 13.8895900677599720 
3  45.8794471255417450 25.2224489680994033 15.8820586780104948 
4  56.5536805700466730 30.8599451390882572 17.7891887400829454 
5  62.0946907412422107 34.4053377263977680 19.6360900260104644 
10 65.5434280496591840 37.0591533853929107 29.9107276797517088 
15 66.3036335404902531 37.5597601427863666 42.5936351248714118 
20 67.1529935901520077 38.1175831389633544 56.3495160940938647 
25  67.9670887224138056 38.6538868920041594 69.2085621375481992 
30  68.6434584860992204 39.1006676479843307 79.6580758797434784 
35 69.1416085841688926 39.4304184865951300 87.2220350799413211 
40 69.4767925880774584 39.6526229458336062 92.2496849929260350 
45 69.6888002622855112 39.7933057889316332 95.4044112768102935 
50 69.8176894863002106 39.8788847025348758 97.3128096721919320 
55 69.8941666570641332 39.9296817602021932 98.4417897148602919 
60 69.9388918847710954 39.9593951068871366 99.1008750069828039 
65 69.9648264987910836 39.9766269819894476 99.4826635148563128 
70 69.9797909566821374 39.9865706090160842 99.7028273461452984 
80 69.9933466703957378 39.9955785990514202 99.9021824887235254 
90 69.9978128089895222 39.9985465136219318 99.9678459183538876 
100 69.9992813426285636 39.9995224187655580 99.9894351820121728 
150 69.9999972506609538 39.9999981729340988 99.999959583014500 
200 69.9999999894831860 39.9999999930110804 99.999999845396316 
250 69.9999999999597690 39.9999999999732622 99.999999999408586 
300 69.9999999999998438 39.9999999999999006 99.999999999997726 
400 70 39.9999999999999930 99.999999999999914 
450 70 39.9999999999999930 99.999999999999914 
600 70 39.9999999999999930 99.999999999999914 
1000 70 39.9999999999999930 99.999999999999914 
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TABLE 6. Runge-Kutta numerical results for the parameters of figure 3a 








x(t) 


y(t) 


z(t) 








CONOTKBRWNF oje 


10 
10.1714712848859499 
10.9124848761110318 
12.8068331298185180 
16.9436588893086012 
25.1301967290895085 
39.7909103204179218 
62.7767725920558064 
92.5591664654720746 
122.198951348776703 
142.847672798411224 
150.446114008380448 
147.244649778317012 
138.139423792913192 
127.376912148014682 
117.510847795193698 
97.9345442737386236 
100.225511128516487 
100.324803966953454 
99.9627565982999756 
99.9951714825362216 
100.005620192871206 
99.9995595698490406 
100.000092728494863 

99.999994637777632 
99.999999997910676 
99.999999999999970 
99.999999999999970 
99.999999999999970 
99.999999999999970 


10 
10.0495808974886156 
10.2618637589493673 
10.7743506720148492 
11.7682289368287966 
13.4320617894971832 
15.9707016143782462 
19.7043297675483480 
25.1648483051793406 
32.9975177648245364 
43.5068883 702286514 
56.0273245279359600 
68.7981196794774804 
79.7804871866895270 
87.7486513562281090 
92.6528941295738804 
96.8955070791911482 
98.8900132755399568 
99.9766952411014672 
99.9830093531780762 
99.9847286077376226 
100.000188625448374 
100.000146478656504 
100.000000220496672 
99.999996752820820 
99.999999999550694 

100 

100 

100 

100 


10 
21.0475841767081917 
41.5883314435775162 
73.6303030796512418 
111.718258774337031 
143.545537317488026 
160.325630375715463 
161.422834286784877 
150.505268930363457 
132.918965238387841 
114.442807944181226 
99.2378765786435224 
88.8557971126201808 
83.0599633098661769 
80.9498704441656970 
81.5358933618375518 
96.6979266202355916 
100.442281878854089 
99.7642768126507634 
99.9406334539491184 
100.012627291336770 
99.9980692966264968 
99.9988396803452418 
100.000000729379607 
100.000003484676597 
100.000000000103512 

99.999999999999986 
99.999999999999986 
99.999999999999986 
99.999999999999986 





4.5. The results of fourth order Runge-Kutta method. 


5. CONCLUSION 


With reference to Tables 1, 2, 3 and 4, and Figures 1, 2, 3 and 4, as 
the growth rate of fish species approaches to their harvesting rate and 
as all other parameters involved in a model became relatively small,the 


population peaks (maximum population before attaining the steady state) 


and the equilibrium population level for each fish species increases. The 
fourth order Runge-Kutta numerical solutions matched with the analytical 


solutions of the model. 
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